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ABSTRACT 

We develop an improved method for tracking the nuclear flame during the deflagration phase of a Type la 
supernova, and apply it to study the variation in outcomes expected from the gravitationally confined detona- 
tion (GCD) paradigm. A simplified 3-stage burning model and a non-static ash state are integrated with an 
artificially thickened advection-diffusion-reaction (ADR) flame front in order to provide an accurate but highly 
efficient representation of the energy release and electron capture in and after the unresolvable flame. We 
demonstrate that both our ADR and energy release methods do not generate significant acoustic noise, as has 
been a problem with previous ADR-based schemes. We proceed to model aspects of the deflagration, particu- 
larly the role of buoyancy of the hot ash, and find that our methods are reasonably well-behaved with respect to 
numerical resolution. We show that if a detonation occurs in material swept up by the material ejected by the 
first rising bubble but gravitationally confined to the white dwarf (WD) surface (the GCD paradigm), the den- 
sity structure of the WD at detonation is systematically correlated with the distance of the deflagration ignition 
point from the center of the star Coupled to a suitably stochastic ignition process, this correlation may provide 
a plausible explanation for the variety of nickel masses seen in Type la Supernovae. 

Subject headings: hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: general 
— white dwarfs 



deflagration stage. This formalism will be used for studying 
a variety of features of all of the above paradigms in future 
work. Nucleosynthesis of species produced as a result of elec- 
tron captures provides a very important observational con- 
straint on supernova models, especially the pure-deflagration 
scenario. For th is reason, and becau se methods are available 
in the literature (Gamezo et al.ll2005!) . we have included elec- 
tron capture and neutrino emission in the energetic treatment 
to capture its effect on the hydrodynamics . Description of this 
method incorporates our previous work dCalder et al.l|200'7b 
detailing the nuclear processing of '^C and '^O by a flame 
front and the evolution of the resulting ash. A method for in- 
tegrating this simplified 3-stage energy release with an artifi- 
cially broadened flame is described in Section|2] The acoustic 
properties of this method are discussed in Section [3] where it 
is shown that the front emits very little acoustic noise. This is 
important for reducing spurious seeding of the strong hydro- 
dynamic instabilities present during the deflagration phase. 

We have chosen in this work to initially pursue simulations 
of the deflagration phase in GCD because it provides a more 
direct demonstration of the buoyancy character o f the flame 
bubble and current work on this mechanism (iPlewa et al.l 



L INTRODUCTION 

It is widely believed that in a type la supernova explosion, 
a WD near the Chandrasekhar limiting mass is disrupted by 
a thermonuclear runaway in its interior, and more precisely 
that a subsonic deflagration must precede any detonation (see 
[Hillebrandt & Niemeyer 2000 and references therein). The 
current leading paradigms for how the deflagration of a WD 
takes place, and how this leads to astrophysical properties that 
match observations, are generally termed (1) pure deflagra- 
tion, (2) deflagration detonation transition (DDT), (3) pulsa- 
tional detonation, and (4) gravitationally confined detonation 
(GCD, Ple wa et al.l l2004y In all but the first of these, a sub- 
sonic deflagration phase expands the WD, lowering its den- 
sity, and a subsequent supersonic detonation then incinerates 
the remainder of the star. Among the remaining three, the 
process that is proposed to ignite the detonation is very differ- 
ent, though it is crucial to determine how much expansion can 
occur prior to the detonation in order to predict the variation 
of nickel mass and therefore brightness among the observed 
Type la's. 

A primary purpose of this work is to set out a numerically 
efficient method for modeling the nuclear energy release in 
the flame front that propagates via heat diffusion during the 
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l2004tlCalder et al.l2004tlPlewal2007tlRB"pke et al.l2006bl) can 

benefit from a concise parameter study. In this mechanism, 
the strong eruption of a rising flame bubble through the sur- 
face creates a wave of material traveling over the surface that 
collides at the point opposite breakout, compressing and heat- 
ing unburnt surface material until detonation conditions are 
reached. We therefore proceed in section |4] to discuss our 
setup for simulating the deflagration of the star, in which our 
principle hypothesis is that the first flame ignition point is rare 
and therefore the deflagration phase is dominated by a single 
flame bubble. Some perspective is given with respect the con- 
ditions expected to be present in the WD core at this time, 
and we describe the progression of the burning in the simula- 
tion, including a survey of the effects of simulation resolution. 
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Finally, in section |5] we present the results of simulations in 
which the ignition point of the flame is placed at various dis- 
tances from the center of the star. We find that the density 
of the star at the time when the GCD mechanism predicts an 
ignition of the detonation, and thus the mass that will be pro- 
cessed to Fe group elements, is well correlated with the offset 
of the initial ignition point. This parameter study also serves 
as a touch-point for future larger-scale simulations of this 
mechanism in three dimensions (Jordan etal. 2007), which 
are essential forjudging its viability (Ropke et al. 2006b). We 
summarize and make some concluding remarks in section|6] 

2. BURNING MODEL FOR A CARBON OXYGEN WHITE DWARF 

There are two fairly different methods of flame-front track- 
ing used in contemporary studies of WD deflagrations. Use 
of a front-tracking method is necessary because the physical 
thickness of the flame front is unresolvable in any full-star 
scale simulation, with the carbon consumption stage being 
10""* to 10^ cm th ick for the density range important in the star 
(ICalder et al.l200 7). The method presented here is based upon 
propagating a reaction progress variable with an advection- 
reaction-diffusion (ADR) equation, and can be thought of as 
an artificially thickened flame, because the real flame is also 
based on reaction-diffusion on much smaller scales. This 
type of method has been used in many previous simula- 
tions of_the_WD deflagration, both in full star simulations 
(e.g. lGamezoTt al. 2003; Cald er et al.ll2004l: |Plewa 2007) and 
to study the effect of the Ray leigh-Taylor (R-T) instability 
on a propagating flame front (iKhokhlovl 1 1 9951 : IZhang et"al] 
120071) . Our flame propagation is based heavily on this work, 
and we have made several refinements to the method that 
we will describe in detail below. The other widely used 
method utilizes the level set technique (ISmilianovski et al] 
Il997t iReinecke et al.lll999l: iRopke et al.ll2003h and performs 
an interface reconstruction in each cell based upon the value 
of a smooth field defined on the grid and propagated with an 
adv ection equation actin g in addition to the hydrodynamics. 
See lRopke et al.1 (1200 6b') andl Schmidt etal.1 (12006.) for recent 
deflagration simulations using this method. 

It should be emphasized that the implementation of the 
flame propagation is far from the only difference between 
these approaches, and there is considerable latitude even 
within one of the front- tracking methods. In addition to the 
front-tracking itself two other issues are important. First, the 
energy release of the nuclear burning must be treated, and this 
is typically done in some simplified way for computational 
efficiency. For example, a prominent difference between the 
method presented here and that commonly used with level-set 
is that we include electron captures in the post-flame material 
within our treatment. The second important additional com- 
ponent is what measure is taken to prevent the breakdown of 
the flame tracking method when R-T, and possibly secondary 
instability in the induced flow, is strong enough to drive flame 
surface perturbations on a sub-grid scale. Both methods fail in 
this limit because the scalar field being used to propagate the 
flame is distorted by advection due to strong turbulence. Gen- 
erally this has been overcome by increasing the flame speed 
enough to polish out grid-scale disorder in th e flow field. This 
can, howeve r, be phrased in term s of simple (IKhokhlovl 1 995h 
or complex dSchmidt et al.ll206 ^ laws intended to mimic the 
enhanced flame surface area produced by unresolved struc- 
ture in the flame. We will leave further discussion to separate 
work, but awareness of this difference is important for com- 
paring results of the two methods. 



2.1. ADR Flame-front Model 

Generally, an ADR scheme characterizes the location of a 
flame front using a reaction progress variable, 0, which in- 
creases monotonically across the front from (fuel) to 1 (ash). 
Evolution of this progress variable is accomplished via an 
advection-diffusion-reaction equation of the form 

^ + v-Vc^=kV^c^+-R{4,), (1) 
at T 

where v*is the local fluid velocity, and the reaction term, R{(t)), 
timescale, r, and the diffusion constant, k, are ch osen so that 
the front propagates at the desired speed. Vl adimirova et alj 
(2006) showed that the step-function reaction rate widely in 
use led to a substantial amount of unwanted acoustic noise. 
They studied a suitable alternative, the Kolmogorov Petrovski 
Piskunov (KPP) reaction term which has an extensive history 
in the study of reaction-diffusion equations. In the KPP model 
the reaction term is given by 

/?(</.)= i0(l-<^). (2) 

The symmetric and low-order character of this function gives 
it very nice numerical pr operties, leading to amaz ingly little 
acoustic noise. Following Vladimirova et an(l2006i) . we adopt 
K, = sbAx/ 16 and r = b Ax /16s, where Ax is the grid spacing, 
s is the desired propagation speed, and b sets the desired front 
width scaled to represent approximately the number of zones. 

The KPP reaction term, however, has two serious draw- 
backs. Formally, the flame speed is only single valued for 
initial conditions that are precisely zero (and stay that way) 
outside the burned region (Xin 2000), which cannot really be 
effected in a hydrodynamics simulation. This can lead to an 
unbounded increase of the propagation speed, which is pre- 
cisely the property we wish to have under good control. Sec- 
ondly, the progress variable cf) takes an infinite amount of time 
to actually reach 1 (complete consumption of fuel). While not 
a fatal flaw like the flame speed problem, this is a problem for 
our simulations in which we would like to have a localized 
flame front so that fully-burned ash can be treated as pure 
NSE material. 

Both of these drawbacks can be ameliorated by a slight 
modification of the reaction term ( Asida et al.. 2007) to 

^(0)={(0-eo)(l + (3) 

where < £[),ei <C 1 and / is an additional factor that de- 
pends on eo and ei and the flame width so that the flame speed 
is preserved with the same constants as for KPP above. This 
"sharpened" KPP (sKPP) has truncated tails in both directions 
(thus being sharpened), making the flame front fully local- 
ized, and is a bi-stable reaction rate and thus gives a unique 
flame speed (iXinii,200 0). The price paid is that since R{(f>) = 
for (j> <Q and (j>> 1, (O is discontinuous, adding some noise 
to the solution. Since the suppression of the tails is stronger 
for higher eo and ei , we adjusted eo and ei so that for a particu- 
lar flame width we could meet our noise goals. The parameter 
values used in the simulations presented in this work were 
eo = ei = 10"^, / = 1 .309, and b = 3.2. The noise properties of 
these choices are discussed in section[3] 

Diffusive flames are known to be subject to a curvature ef- 
fect that affects the flame speed when the radius of curva- 
ture is similar to the flame thickness, a frequent circumstance 
with modestly-resolved flame front structure. In testing, the 
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curvature effect of the step-function reaction rate proved sur- 
prisingly strong, likely due to the exponential "n ose" that the 
flame front possesses (l yiadimirova etai] I2OO6'). Both KPP 
and our sKPP show significantly better curvature properties. 
Due to the necessary discussion of background and the size 
of the study supporting this conclusion, this t opic will be dis- 
cussed in detail separately dAsida et alj|2007h . 

2.2. Brief Review of Carbon Flame Nuclear Burning in 
White Dwarfs 

In previous work dCalder et al.ll200'7l) . we performed a de- 
tailed study of the processing that occurs in the nuclear flame 
front and th e ashes it leaves . It wa s shown that, as discussed 
previously (lKhokhlovlll983L |1991|) . the nuclear burning pro- 
ceeds in roughly 3 stages: consumption of '^C is followed by 
consumption of '^O on a slower timescale, which is in turn 
followed by conversion of the resulting Si group nuclides to 
Fe group. Most of the energy release takes place in the '^C 
and '^O consumption steps, and at high densities the result- 
ing material contains a significant fraction of light nuclei (a, 
p, n) and is in an active equilibrium in which continuously 
occurring captures of the light nuclei are balanced by their 
creation via photodisintegration. Initially, the heavy nuclei 
are predominantly Si group, this is termed nuclear statistical 
quasi-equilibrium (NSQE), which upon conversion of these to 
Fe group becomes nuclear statistical equilibrium (NSE). Each 
of these states is reached on a progressively longer timescale, 
and the importance of distinguishing the last lies in the dis- 
parity in electron capture rates between Si and Fe group based 
equilibria. The energy released by our scheme at a given den- 
sity has been dire ctly verified within a few percent against 
those tabulated in ICalder et al.l ( l2007h and against an addi- 
tional direct NSE solution. 

O ur methods build h eavily on t hose of|Gamezo et al.l (l2005h 
and iKhokhlovl (flMh (see also iKhokhlovl l2000r which is 
used t hroughout their fam i ly of recent deflagration calcula- 
tions ( Gamezo etaT]l2003l |2004 |2005). The principal dif- 
ferences, other than the use of the sKPP reaction term, are 
that we use the predicted binding energy, qf (see definitions 
below), of the final NSE state rather than using qa^eip, T, Ye) 
with the current density and temperature, p and T, and we sep- 
arately track '^C and "'O consumption. These are described 
in detail below. Finally, the method presented here is entirely 
different from that used by Plewa (2007), which effectively 
"freezes" the NSE at the state produced in the flame front, ne- 
glecting the additional energy release as the light nuclei are 
recaptured. 

2.3. A Quiet Three Stage, Reactive Products Flame Front 

In incompressible simulations, the progress variable in an 
ADR front tracking scheme is typically used to parameterize 
the density or density decrement. An analog in compressible 
simulations is to release energy in proportion to (/>. This sim- 
ple idea becomes somewhat complicated in a situation like the 
WD, where the burning (and therefore energy release) occurs 
in multiple stages whose progress time scales vary by orders 
of magnitude during the simulation. A further complication 
is created by the dynamic NSE state of the ash, such that the 
energy release depends on the physical conditions (density) 
under which the flame is evolving. 

The ethos we have implemented here is that processes that 
occur on scales that are unresolved by the artificially thick- 
ened flame should have their energy release counted towards 



the overall energy that is smoothly released by the progres- 
sion of the flame variable 0. This approach is accomplished 
by defining additional progress variables that follow the ADR 
variable cj) and that govern the energy release. Such a com- 
plex scheme is necessary for the nuclear burning in the WD 
because, as shown in Calder et al.l (12007 ) . the conversion of Si 
to Fe group that occurs over centimeters near the core, occurs 
over kilometers in the o uter portion of the star In previous 
work (ICalder et al.l2007h . we presented a method for integrat- 
ing energy release with an ADR flame. Those prescriptions 
were an early version of what is presented here, and are su- 
perseded by the method presented below. Note in particular 
that the functional meaning of 02 and 03 have changed some- 
what because the ash state is able to evolve regardless of the 
value of the progress variables. Also, the use of surrogate nu- 
clei described in that work has been abandoned in favor of the 
direct use of scalars described below. 

We define three progress variables, which represent irre- 
versible processes. These three variables start at in the un- 
burned fuel and progress toward 1, representing 

01 Carbon consumption, conversion of C to Si group 

02 Oxygen consumption, conversion of O to Si group 

03 Conversion of Si group to Fe group. 

We keep strictly 0i > 02 > 03, but all three are allowed to 
have values other than zero and 1 in the same cell. It is most 
useful to think of material in a cell as being made up of mass 
fractions of 1 -0i of unburned fuel, 0i -02 of partially burned 
(no carbon) fuel, 02 of NSQE material of whic h 03 has had 
its Si group elements consumed. As shown by ICalder et al.l 
given sufficient resolution, all these stages are, in fact, 
discernible as fairly well separated transitions. However, with 
an artificial flame, a real transition from fuel to final ash that 
occurs in less than one grid spacing must be spread out over 
several. 

We now describe how these auxiliary progress variables 
track the flame progress. In our case the evolution of 0i is 
set directly by the artificial flame formalism described above, 
01 = 0. Thus the noise properties of the artificial flame it- 
self are inherited by the energy release scheme. The connec- 
tion between the energy release and the ADR flame tracking 
comes entirely through this equality, and so coupling the fol- 
lowing energy release methodology to other available front- 
tracking methods appears quite practicable. 

We evolve a number of scalars which, in the absence of 
sources, satisfy a continuity equation, 

^ + V-(Qpy) = 0, (4) 

where Q is the scalar under consideration. The flame vari- 
able above is one such scalar, and our additional progress 
variables are also treated as such. The other scalars we utilize 
directly represent physical properties of the flow; they are the 
number of electrons, the number of ions and nuclear binding 
energy per unit mass or baryon, respectively: 

i'-Ef:^" (5) 

^"=E^^" (7) 
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where / runs over all nuclides, Z,- is the nuclear charge. A, 
is the atomic mass number (number of baryons), and Ei,j = 
(Zinip - Nim„ - mi)c^ is the nuclear binding energy where Z,, 
nip, Ni, and m„ are the number and rest mass of protons and 
neutrons respectively, so that positive is more bound. The 
mass fractions Xj are not treated in our simulation, and are 
used here only to define these properties, though ©-(iTll sat- 
isfy (|4|l by virtue of being linear combinations of the mass 
fractions, which themselves satisfy (|4]i. Defining our flame 
model is then a matter of setting out the source terms for the 
various scalars (pi, (pi, </>3, Y^, Fion, and q, and relating these to 
the energy release. 

For a given (pi and (p2 we define, for notational conve- 
nience. 



Xc = {l-cpi)Xl 

Xo = (l-02)(1- 
XM^ = {<pl-<p2)^l, 



(8) 
(9) 
(10) 



where X^ is the initial carbon fraction. These can serve as 
approximate abundances, though the real abundances in these 
stages have several additional important species. Since this 
notation can be misleading, we again emphasize that abun- 
dances are not being tracked in our simulation, the material 
properties used in the EOS are derived directly from Yg and 
Fion, discussed further below. The remaining mass fraction 
of material 02 is considered to be in NSQE or NSE, so that 
02+^c+-'^o+-'^Mg = 1- We define this "ash" material to have 
binding energy ^ash and electron fraction Fg ash and ion number 
i'ion.ash such that, 



q = 02?ash +Xcqc+ Xoqo + XMgqUg 



(11) 



and similarly for the other quantities. To again clarify our 
notation, ^{c,o,Me} are the actual binding energies of '^C, ""O 
and ^"^Mg, being used here to approximate the binding energy 
of the intermediate ash state. The final scalar, 03, represents 
the degree to which the ash has completed the transition from 
Si-group to Fe-group heavy nuclei, and is used to scale the 
neutronization rate as described below. Thus after material 
has expanded and is no longer a-rich, Xsi-group ~ 02 - 03- 

From the quantities q and the local internal energy per 
mass, £, it is possible to predict the final burned state if 
density, p, or pressure, P, were held fixed for infinite time 
and weak interactions (e.g. electron captures) were forbid- 
den (constant Ye); this gives the NSE state. Our equations 
and formalism for NSE, w hich include plasma Coulomb cor- 
rections, were described by lCalder et al.l (l2007l) . Using these, 
the abundances and therefore average binding energy of the 
NSE state can be found for a given p, T, and Ye, resulting 
in qNSEip,T,Ye). Fion.NSE, as well as the Coulomb coupling 
parameter F = Z^^^e^(4-TTne/3y^^/kT, where e is the electron 
charge, n^ is the electron number density and k is Boltzmann's 
constant, are similarly determi ned. For i nternal energy, £, we 
follow the convention of Tim mes & Swes ty (2000), which ex- 
cludes the rest mass energy of the (matter) electrons. 

The final burned state at a given p and Y^ can be found by 
solving 

£-q = £(Tf)-qj,sE(Tf) (12) 

for Tf. The number of protons, neutrons and electrons is the 
same in both states, so that the rest masses cancel in this equa- 
tion. We will denote qj- = ^nse(7/)> as the solution to this or 
the corresponding isobaric equation below. For the sake of 
computational efficiency, this solution is accomplished via a 



table lookup in a tabulation of q/ip, Ye,£- q). Since the flame 
is quite subsonic, it is also useful to be able to predict the 
NSE final state for the local P. This can be accomplished by 
solving, at a particular P and Y^, 



_ P _ P 

= f(7»-^NSE(7» + — — , 
P P(Tf) 



or, more naturally. 



n-q = n{Tf)-q(Tf), 



(13) 



(14) 



where H is the enthalpy per unit mass. This leads to a sim- 
ilar tabulation of qf{P,Yi,,T-l- q). We denote the solution of 
( fT2] i as isochoric and that of (fT4l i as isobaric. While the iso- 
baric prediction of the final state must be used within the 
thickened flame front, away from the flame front we would 
like to use the isochoric result to avoid undue interference 
with the hydrodynamic evolution. This necessitates a hand- 
off when the material is nearly fully burned. We wish the 
hand-off to occur at a high enough that the difference be- 
tween the isochoric and isobaric results are minimized, but we 
introduce a small region where the results are explicitly aver- 
aged in order to minimize the noise generated by the hand-off. 
Thus, where 02 > 0.9999 we use the isochoric estimate, for 
0.99 < 02 < 0.9999 we use a linear admixture of isochoric 
and isobaric, and for 02 < 0.99 we use the isobaric estimate. 
From noise tests and behavior in simulations, these values ap- 
pear sufficient for the current purposes. 

We now have in hand an estimate of the NSE final state 
qf, and its temperature Tj. Calder et al. (2007) evaluated the 
timescales for progression of the burning stages from self- 
heating calculations, as functions of Tj. The progress vari- 
ables are then evolved according to'' 



(15) 
(16) 



'nvSQE(7}) 

02-03 



0-i = - 



and the binding energy of the ash material according to 

?ash = ■i-qf+ 7^ ■ (17) 

02 tnsqe(-'/) 

This is in fact implemented conservatively in the finite differ- 
ence form 



«+i 

^ash 



?:sh+^^^AO+^/02Af 
tnsqe(^/) 



/(02 + 02Af) , 

(18) 

where At is the timestep. The ion number is treated similarly, 
according to 



y _ , yion./-yion 

-'ion. ash ~ , -'ion./"'" . 

02 tmsqeUj) 



(19) 



with a similar finite differencing. Neutronization (mainly 
electron captures) is implemented by applying 

t,.sh = <p3Yeip,Tf,Ye) ■ (20) 

Our calculation of Y^ is described in ICalder et al.l (l2007h and 
utilizes 443 nuclides in the NSE calculation including all 

' Here X denotes a Lagrangian time derivative, tiiougii since our code is 
operator-split between tiie liydrodynamics and source terms, tlie implemen- 
tation is a simple time difference. 
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available rates from lLanganke & Martmez-Pinedol(l2000h . Fi- 
nally q is recalculated with the new values of (j>i, (f>2 and ^ash 
using ( fTTT ) and the energy release rate is 

• 2 

Enuc = -^-hlYeNAC {mp + Mg- Mn) + e^] ■ (21) 

Here ei, is the energy loss rate from radiated neutrinos and 
antineutrinos, and is calculated similarly to Y^. 

Our description here has been fairly algorithmic, for the 
sake of clarity of implementation. It is possible, however, to 
use eq. ([TsTl-lfTTli along with the definitions (ISTl-dTOll and ( fTTT l 
and ( |2TI ) to obtain straightforward Lagrangian source terms. 

3. QUANTIFYING ACOUSTIC NOISE FROM THE NUMERICAL 
FLAME FRONT 

The noise generated by the model flame may influence the 
outcome of a deflagration simulation by seeding spurious fluid 
instabilities. Quantifying noise, determining the sources of 
noise, and minimizing noise are therefore necessary steps in 
the development of a robust flame model. To this end, we 
pe rformed a suite of simple test simulations following those 
of IVladimirova et all (l2006l) . The simulations presented here 
are for the sKPP flame model with eo = f-i = 10"^, the high- 
est values for which the RMS deviation in velocity (see be- 
low) was a few X 10"'* or lower for the density range 10^- 
lO' g cm"-'. We note that simulations of model flames utiliz- 
ing the "top hat" reaction produced considerably more noise, 
~ 0.1 or more RMS velocity deviation. 

3.1. Details of Test Simulations 

The simulations consisted of one-dimensional flames prop- 
agating through 40 km of uniform density material composed 
of 50% and 50% '^O by mass. The simulations had 
a reflecting boundary condition on the left side and a zero- 
gradient outflow boundary condition on the right with the 
flame propagating from left to right. The flame was ignited by 
setting the left-most 5% of the domain to conditions expected 
for fully burned material, with the transition to unburned fuel 
described by the expected flame profile. This method of igni- 
tion approximated what would have resulted from letting the 
flame burn across the ignition region. 

The choice of boundary conditions allowed material to flow 
to the right and off the grid as the flame propagated and more 
of the domain consisted of (expanded) ash. The densities, 
sound speeds, and sound crossing times of the simulation do- 
main are given in Table[T] The simulations were performed on 
domains of 256, 512, and 1024 zones, corresponding to reso- 
lutions of 15625, 7812, and 3906 cm respectiv ely. The sim- 
ulations were performed with the FLASH code (iFryxell et al.l 
I2OOOI: ICalder et al. 2002). which explicitly evolves the equa- 
tions of hydrodynamics. In an explicit method such as this, 
the time step of a given simulation is limited by the sound 
crossing time of the zones. For these simulations, the maxi- 
mum time step was determined by a CFL limit of 0.8, mean- 
ing that the time step allowed a sound wave to cross only eight 
tenths of the zone with the highest sound speed. 

Acoustic noise may be quantified by considering the mag- 
nitude of variations in quantities like pressure and velocity. 
We define the "RMS deviation" of a quantity x as 

devRMs(^) = y^(x2)-(x)V (x) , (22) 

where the averages are taken in space at a given time. In each 
simulation we calculated the RMS deviation of pressure and 




Time Step 

Fig. 1 . — Plot of RMS deviation in pressure and velocity for simulations of 
propagating model flames at a density of lO' g cm"'' performed on simulation 
domains of 256, 512, and 1024 zones. Shown are the RMS deviations as 
functions of the number of time steps, with solid curves indicating pressure 
and dashed curves indicating velocity. In these, the flame speeds were 3.89 X 
10^ cm s"' , the expected flame speed for material of this density. 

velocity in both the fuel and the ash, with fuel defined as the 
region of the domain with 0i < 0.0001 and ash defined as the 
region with 0i > 0.9999. The RMS deviations presented in 
the figures here are for the fuel. 

The densities we considered ranged from w 10^ to 2 x 
10' g cm"-'. Nuclear flame speeds vary extremely with den- 
sity, from « 5 X 1 0^ to 8 x 10^ cm s"' for laminar flame speeds 
at the se densities jTimmes & Wooslevl ll992l;IChamulak et al.l 
120071) . Because of the disparity of flame speeds, we per- 
formed some simulations with a fixed flame speed of 6 x 
10^ cm s"', allowing flames to propagate across the simula- 
tion domains in similar elapsed time. Note the time for the 
flame to cross the domain is shorter than the domain size over 
the flame speed due to expansion across the flame, which it- 
self depends on density. 

Figure [T] shows the RMS deviation in pressure and veloc- 
ity for three simulations at a density of lO' g cm"^ plotted 
as functions of the time step number for each simulation. 
Plotting the RMS deviations against the time step number 
is equivalent to scaling the evolution time of the simulations 
by the fraction of the sound crossing time of each simula- 
tion and eliminates the differences in time step size due to the 
different resolutions. In this case the input flame speed was 
3.89 X 10^ cm s"'. 

Figure |2] presents the RMS velocity and pressure devia- 
tions from fixed flame speed simulations performed at p = 
10^,10^,10' gem"-'. The panels present a set of simula- 
tions performed on a uniform simulation domain of 256 zones 
(top), 512 zones (middle), and 1024 zones (bottom). In this 
figure, the RMS deviations are plotted against the simulation 
time. 

3.2. Sources of acoustic noise 

The simulations performed for this study indicate that there 
are two principal sources of noise, transient noise resulting 
from the initial conditions and steady rhythmic noise pro- 
duced as the model flame propagates across the simulation 
domain. We observed noise from these sources in the simula- 
tion results in three forms, described below, one from the ini- 
tial transient and two from the propagation noise. Though not 
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Fig. 2. — Plot of RMS deviations in pressure (solid curves) and ve- 
locity (dashed curves) for simulations of propagating model flames per- 
formed at densities of lO' 10^, and lO' g cm"^^ with a fixed flame speed of 
6 X 10* cm s"' . The simulations were performed on a uniform domain with 
256 zones (top panel), 512 zones (middle panel), and 1024 (bottom panel) 
zones. 



TABLE 1 
Simulation Properties 



density 


Cs (fuel) 


^soundcrossing 


10' g cm-3 


lO** cm s"' 


lO-'s 


2 


9.04 


4.42 


1 


8.06 


4.96 


0.5 


7.17 


5.58 


0.3 


6.58 


6.08 


0.1 


5.50 


7.27 


0.05 


4.82 


8.30 


0.03 


4.40 


9.09 


0.01 


3.59 


11.1 



shown here in detail, by comparing these metrics for the mul- 
tistage flame to a single-stage flame with comparable static 
energy release, we confirmed that the multistage scheme adds 
no significant noise. 

The most obvious feature of the simulations is a large 
amount of noise early in the simulations. This may be ob- 
served in all of the RMS deviation figures as the large devi- 
ations on the left-hand side (early time region) of the plots. 
This transient results from the method of initiating the burn- 



ing. The burned region is created by setting Fion./, and qj to 
the appropriate values for the given density and fuel compo- 
sition and the profile of the flame to 

0iW = i [1.0-tanh [(x-xo)/L]] (23) 

where L = bAx/2 (IVladimirova et alJl2006l) . 02 = </>!, ^3 = 0, 
and here xq is the initial position of the flame front, 5% of the 
distance across the domain. This prescription produces initial 
conditions that depart slightly from the relaxed result obtained 
after some evolution, and the relaxation occurring during the 
initial evolution produces the large amounts of noise observed 
early on. The perturbation in this case is a large sound pulse 
that propagates across the domain. 

Observation of the RMS deviation curves for a particular 
density in Figure |2] indicates that the duration of the transient 
noise depends on resolution of the simulation. This occurs be- 
cause the width of the thickened flame is set by the resolution 
of the simulation grid, and the duration of this pulse is set by 
the flame self-crossing time. The wider flames at lower res- 
olutions produce an initial transient sound wave that has the 
correspondingly longer wavelength thereby taking the corre- 
spondingly longer time to all propagate out of the domain. As 
an example we consider simulations of Figure [T] performed 
with p = 10"^ g cm"-' and s = 3.89 x 10^ cm s"' . The observed 
times for the transient pulse to pass completely across the grid 
were 0.016, 0.026, and 0.047 s for resolutions of 1024, 512, 
and 256 zones, respectively. These times were measured by 
observing the pressure wave propagate across the simulation 
domain and agree well with the duration of the transients ob- 
served in the deviations. The times are very consistent with 
the flame self-crossing time, 4-Ax/s, which is also the time for 
the flame profile to come into equilibrium, and therefore for 
the burning rate to stabilize. Thus the duration in number of 
time steps, 4Ax/s/(Q.^Ax/cs) ~ 1000, is similar for all three 
resolutions, as seen in Figure [T] 

As the flame profile moves across the regular underlying 
grid, the slight changes in the profile due to the spatial quanti- 
zation lead to a small, rhythmic variation of the burning rate. 
This produces a pressure wave train propagating out through 
the fuel with a specific form characterized by the quantization 
and with a period determined by time for the flame profile to 
shift by one zone. This leads to the high frequency oscilla- 
tions readily observed in Figure[T]after the initial transient. In 
Figure |2] this noise may be seen as the small high-frequency 
(period < 0.005 s) features on the curves and is most obvi- 
ous in the p= 10^ g cm"-' curves of the 256 zone simulation 
(top panel). As an example, we consider the simulation at 
p = 10^ g cm"-' with a flame speed of s = 3.89 x 10^ cm s"' 
at our highest resolution, shown in Figure [T] The wave 
propagating through the fuel has an average wavelength of 
6.76 X 10^ cm, which with a sound speed of 8 x 10^ cm s"' 
gives a period of 8.45 x 10""* s. The flame front propagates 
through 16 computational zones in an average of 0.0135 s, 
giving 8.43 x 10"'*s to burn each zone. With the fuel at rest, 
the flame front should move within the computational domain 
at ipfuei/Pash (this is higher than the flame speed due to ex- 
pansion), which gives a zone crossing ti me of 8.37 x 10"^ s , 
using the expansion for this density from lCalder et al.l (l2007h . 
All these checks prove completely consistent. 

The regular oscillating pattern (with period 0.1s and larger) 
of the noise visible in Figure |2] originates from the variation 
in size of the region of the emitted wave train over which 
the deviation is being averaged. From the discussion above 
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this wave train has a wavelength given by equating the sound 
crossing time of the disturbance with the time for the flame 
to cross a single zone, A/q = Ax/vy. Since the sound field 
is otherwise essentially flat, averaging over an integer number 
of these wavelengths will give about the same result. Thus 
we expect a regular pattern in the noise measured at a period 
determined by the time for the averaging interval to shrink by 
one wavelength. The interval is shrinking at the same speed 
that the flame is moving across the computational domain, 
given above, and dividing the wavelength of the emitted train 
by this gives a period of P = A/v/ = CsAji:/(ipfuei//5ash)^- This 
relation reproduces the linear dependence on resolution seen 
in the results, and shows that the dependence on density enters 
through both the sound speed and the expansion factor We 
have confirmed that this relation reproduces the periods in the 
figures, e.g. P = 0.24 s for s = 6 x 10^ cm s"\ p= 10'' g cm"-' 
and our coarsest resolution. This result represents two bumps 
in the noise figures because we are taking the RMS deviation, 
losing the sign. 

Finally, we note that as the flames approached the edge of 
the simulation domain and most of the fuel on the domain 
had burned, the magnitude of both the high frequency oscil- 
lations and the regular pattern in the deviation increased (as 
may be observed on the right hand side of the curves, espe- 
cially at the lowest density, 10^ g cm"-'). This increase occurs 
because what little fuel remains samples the region very near 
the flame, which is expected to have the most noise. 

4. THE PROGRESSION OF A FLAME FROM SINGLE-POINT 
INITIATION 

It is useful to describe with some detail the progress of 
events involved in the GCD mechanism, and how our sim- 
ulation setup captures these events. In the centuries before 
the la event, when the WD has accreted enough matter to ig- 
nite carbon burning in the center, there is an expanding con- 
vective carbon-burning core (see e.g. Wooslev et al. 2004). 
This state is already a runaway, because the temperature at 
the center will continue to monotonically rise. Ignition oc- 
curs during this convective phase when the local heating time 
Theat — CpT /e„uc, whcrc cp is the specific heat at constant pres- 
sure and enuc is the nuclear energy deposition rate, becomes 
shorter than the eddy turnover time Tgdd — 10-100 seconds. 
At this point the burning runs away locally, the '^C and '^O 
fuel converts entirely to Fe-peak elements and a flamelet is 
born. 

While the rate of formation of ignition points is unclear, 
it is believed that ignition of local fla melets in the core of 
the WD is a fairly stoch astic process (IWooslev et al.l l2004t 
IWunsch & Wooslevll2004l) . For this study we will work under 
the hypothesis that the ignition conditions are rare at the time 
the ignition occurs. This means that the ignition grows from 
a small (< 1 km) region somewhere in the first temperature 
scale height (~ 400 km) near the center of the star, and the 
second igni tion is long enough a fter this (> 1 sec) to be unim- 
portant (see lWooslev et al.]|2004l for a discussion of how these 
scales arise.) Thi s picture is representati ve of the ignition con- 
ditions found bv lHoflich & Stein (2002) in their study of the 
pre-runaway phase, but is somewhat in contrast to the conclu - 
sions of some of the above work (e.g. IWooslev et al.ll2004h . 
and is essentially t he opposite hypothesis to that taken by 
iRopke et al.l(l2006ah . Single-point ignition is plausible within 
the current uncertainties and is quite useful for our purpose 
here of understanding the dynamics of a flame bubble and 
characterizing our numerical methods. 



4.1. General Simulation Setup 

In order to simplify our study of single-bubble dynamics we 
begin our simulation with no velocity field in the stellar core. 
This is, in fact, a poor approximation to reality because the 
typical outer scale ve locity in the convec tive core is expected 
to be - 100 km s"' dKuhlen et_al]|2006|), which IS compara- 
ble to the laminar flam e propag ation speed in this part of the 
star (Timm es & Woosleylll992f IChamulaket al. 2007). This 
means that the strongest initial source of perturbations on the 
flame surface, and therefore seeds of the later R-T modes, is 
likely to be the turbulence in the convective core flow field. 
The convection field is, however, not strong enough to de- 
stroy the flamelet once it is born. We feel neglection of the 
the convective flow field in this initial study justifiable for two 
reasons. First, we would like to understand the dynamics of 
the bubbles and flame surface near the core first without the 
additional complication of the turbulent flow. Second, it will 
be challenging to interpret the effects that a turbulent field in 
two dimensions subject to the imposed axisymmetry might 
produce. Any off-axis feature acts effectively as a ring, an ef- 
fect that cau ses enough difficul ty even with static initial con- 
ditions. Also lLivne et al.l (l2005|) find that the general outcome 
of off-center ignitions is not strongly effected by the presence 
of a convection field. 

We perform two-dimensional axisymmetric simulations 
with the FLASH adaptive-me s h hy drodynamics code 
jFrvxeU et al.1 l2000t ICalder et al.l l2002h . We be gin our 
simulations with at 1.38 Mq WD with a uniform composition 
of equal parts by mass of '^C and '^O. This model has a 
central density of 2.2 x lO** g cm"-', a uniform temperature 
of 4 X 10^ K, and a radius of approximately 2,000 km. A 
spherical region on the symmetry axis is converted to burned 
material with (pi given by eq. (|23] | with x = \r- Foff | and 
xo = ?"bub, 4>2 = <l>i, and 03 = 0, where Foff is the location of the 
center of the ignition point. The density is chosen to maintain 
pressure equilibrium with the surrounding material. Thus the 
radius of the flame bubble, rbub, is the approximate location 
of the (pi = 0.5 isosurface, and all simulations in this paper 
begin with a spherical bubble of radius 16 km. This is the 
smallest bubble that is reasonably well resolved (having (pi 
very close to 1 at the center) at 4 km resolution, that used 
in the parameter study presented in section |5] The basic 
parameters and some results are listed in Tabled and will be 
discussed below. 

Our adaptive mesh refinement has been chosen to capture 
the relevant physical features of the burning and flow at rea- 
sonable computational expense. We choose to refine on strong 
density gradients everywhere, and strong velocity gradients in 
the burned material. At the beginning of the simulation all of 
the star is resolved to 16 km resolution regardless of the max- 
imum allowe d resolution for reasons of hydrostatic stability 
(lPlewall2007h . In regions with p<5xlOgcm'' refine- 
ment is not requested, which includes all of the region out- 
side the star Regions where the flame is actively propagating 
(0.1 < (p < 0.9, and non-trivial flame speed) are required to 
be fully refined in order to properly propagate the flame front. 
In the interest of limiting computational expense, the refine- 
ment is limited to a finest resolution of 32 km outside a radius 
of 2500 km, which is above the the surface of the star in all 
cases treated here. This constraint will be relaxed in future 
work, but the bulk kinetic motion of the surface flow is not 
expected to be affected by this choice. 
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Fig. 3. — Critical wavelength for flame surface perturbations Ar = 67rs^/Ag 
in the initial WD. 

4.2. Stages of Flamelet Evolution 

The major stages in the evolution of the bubble can be 
roughly understood by comparing the bubble's size to the 
critical wavelength for the unstable Rayleigh-Taylor growth 
of flame surface perturbations. We define Af = 6TTs'^/Ag 
(iKhokhlovl 1 9951) . where g is the local gravity, s is the laminar 
flame speed, and A = (/Ofuei-pash)/(/Ofuei + /Oash) is the Atwood 
number Below A^, s is high enough that perturbations in the 
flame surface can be "polished out" by burning, so that the 
surface remains laminar and simple. Above A^., R-T is strong 
enough that perturbations will grow instead, leading to a com- 
plex, disordered^ame structure down to a scale of order A^. 
?See 'Khok hlovlfT995l for an extensive discussion.) As shown 
in Figure [3] Ac drops quickly with radius, mainly due to the 
increasing gravity as one moves away from the center of the 
star, and after this due to the falling flame speed. We take 
'"bub < Ac in all cases, so that our assumpti on of a sphe rical 
bubble at rest is physically justifiable (F isher et al.ll2007t) . 

Initially bubbles can be thought of as moving from lower 
left to upper right ( growing and rising) in Figure |3] (see 
IZingale & Dursil 12007 ). We distinguish three phases of the 
bubble evolution and rise in terms of Ac. Each of these phases 
can be seen in Figure H] where we show the evolution of 
the flame surface with time for two different resolutions: our 
typical resolution, 4 km, and our highest resolution, 0.5 km. 
Initially, while small and near the center with rbub < Ac, the 
flame bubble will grow according to the laminar flame speed, 
roughly keeping a spherical shape. This structure is evident 
at f = 0.2 seconds, where the bubble has akeady grown to a 
radius of about 30 km, twice its original size. 

Eventually as the bubble rises and grows, it will reach the 
point when rbub ^ Ac. For our case this occurs at a little over 
100 km radius. As seen at f = 0.4 sec, the bubble forms a 
R-T roll at approximately its full dimensions, the first scale 
that is unstable. This feature shows some differences with 
resolution, but generally the same kind of feature (a roll) has 
appeared in both simulations. Crossing Ac, the bubble is now 
a R-T unstable flamelet. The subsequent evolution is resolved 
in our simulation until Ac falls below the grid scale. Thus 
we will term the second stage of evolution as the resolved 
R-T stage. The R-T structure becomes unresolved at approx- 
imately 350 and 500 km radius for 4 and 0.5 km resolution 
respectively. Thus by t=0.6 seconds, the 4 km simulation has 
already entered the unresolved regime, while the 0.5 km sim- 
ulation has nearly entered it. This is evidenced by the sepa- 
ration between consecutive contours in the progress variable 
caused by strong advection of material within the flame front. 



The rest of the bubble rise and burn is dominated by un- 
resolved R-T burning. The bubble continues to grow, both 
due to burning and expansion of the material under decreas- 
ing pressure. Its top generally reaches the surface at approx- 
imately 1.0 seconds ("breakout"), after which it mostly exits 
the interior of the star and expands strongly to create the flow 
around the stellar surface. A notable difference between our 
evolution and that seen by (iRopke et al.ll2006"bl) is the pres- 
ence of the "stem" below the risin g bubble. This stern is ab- 
sent in non-reactive rising bubbles f Vladimirovall2007h . but it 
is mysterious that it is completely absent in the level-set reac- 
tive flow simulations. This might be related to the fact that the 
initial bubbles used by iRopke et al.l (l2006 bl) are too large and 
far off-center to capture either the laminar or resolved R-T 
evolution stages, but more investigation is necessary. Though 
both simulations in Figure |4] are unresolved at t = 0.8 sec- 
onds, comparing the two provides a good demonstration of 
the enhancement in flame surface facilitated by smaller scale 
structure. 

As mentioned in section |2] any flame tracking method that 
depends on the advection of a scalar field faces difficulties 
in the highly turbulent flow produced in the unresolved R-T 
stage. Measures must be taken to compensate for the finite 
resolution of the simulation, otherwise the flame ceases to 
propagate correctly because the turbulence destroys the nec- 
essary structure for self-propagation of the scalar field. This 
problem is particularly relevant to the simulation of WD de- 
flagratio ns, as has been extensively discussed by previous 
authors (lKhokhlovlll9"9l iGamezo et alll2003t ISchmddt et all 
2006). R-T can create turbulence down to the scale Ac, be- 
low which the flame i s able to p olish out the mixing pertur- 
bations. It is thought dKhokhlovt 19951) that the details of this 
small-scale behavior need not be followed explicitly. If the 
overall dynamics of the flame surface are determined by the 
large scale perturbations, simulations of moderate resolution 
(which may however be only just possible today in 3-d) can 
determine the evolution of the burning front. This is called 
self-regulation, and the open question is: What resolution is 
"enough"? This is currently under debate, and we will argue 
here that we have enough resolution for some purposes and 
can make statements about others based on trends that we see 
w ith resolution. 

ISchmidt et al.l (l2006h attempted to compensate for the 
shortcoming of the limited resolution by artificially enhancing 
the propagation speed of the flame front based on a measure of 
local turbulence. Such complex measures should not be nec- 
essary if self-regulation holds, thus we have taken what we 
consider to be a conservative approach, making the smallest 
possible adjustment to the front propagation speed and evalu- 
ating the inaccuracy directly via resolution study. We enforce 
a minimum flame speed (which therefore acts effectively as an 
enhancement) that is intended to ensure tr-t ^ Tflame where, 
for flame width 5, tr-t ~ y^S/Ag is the characteristic R-T rise 
time and Tflame ~ is the flame self-crossing time. Since our 
flame width is proportional to the resolution, this leads to a 
limit s > a^AgmjAx, where Ax is the resolution, a is a ge- 
ometrical factor that is different in two and three dimensions, 
and m j is a calibrated constant. 

We have determined nif, the constant that determines 
when the flame model "breaks", empirically by running 
two-dimensional fla me-in-a-box simulat ions like those of 
IKhokhlovl (119951) and lZhang etalj (l2007h . We evaluated the 
range of s over which self-regulation holds, in which the burn- 
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Fig. 4. — Stages of bubble growth at different resolutions. The initial radius of the bubble was 16 km centered at 40 km from the center of the star, contours 
are shown for the progress variable <j) = 0.l (green), 0.5 (red), and 0.9 (blue). 



ing rate, mbumed = pL^'''^^ s^ff is determined by the box size, 
L, such that igff = cty/AgL, independent of s. In two dimen- 
sions a = 0.28 and in three dimensions, as found by previ- 
ous authors, a = 0.5. We found that the self-regulated regime 
is bounded above by s ~ \JAgL, corresponding to the re- 
quirement Ac < L, and below by s ~ a-\/A^0.04Ax. In our 
WD simulations we have used the value for a from three- 
dimensional simulations, although we are working in two- 
dimensional cylindrical geometry, and nxf = 0.04. These val- 
ues have been confirmed by preliminary three-dimensional 
flame-in-a-box calculations, which will be discussed in a sep- 
arate work. Using this type of floor on s requires that we 
explicitly turn off the flame at low density. This is done 
smoothly between p= 10^ g cm"^ and 5 x 10^ g cm"^, so that 
s = for densities below this. 



4.3. Resolution Study 



Some properties of the off-center deflagration model that 
we are trying to deduce from our simulations show depen- 
dence on the simulation resolution, while others do not. We 
would like to make statements as much as possible based on 
features that are not influenced by resolution, and where we 
cannot avoid it, account for the dependence in other conclu- 
sions that we draw. Problems with resolution-dependence is 
not entirely unexpected, since, as just discussed, a significant 
amount of our simulation is a priori known to be unresolved. 
In summary, we find that the conditions at the possible deto- 
nation point are fairly insensitive to resolution, for the reso- 
lutions considered, but that the state of the interior of the star 
at a given time during the runaway may only be calculated 
by higher resolution simulations than the 4 km at which our 
parameter study was performed. 

Of foremost interest is the robustness of the gravitationally 
confined detonation (GCD) mechanism, particularly the prop- 
erties in the collision and compression regions opposite the 
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Fig. 5. — Maximum temperature Tmax in tlie lower hemispliere and the 
density at the same point. Shown are results from simulations where the de- 
flagration is ignited by a bubble of flame with a radius of 16km placed 40 km 
off-center and for which the resolution is varied to 8, 4, 2, and 1 km. Ma- 
terial flowing over the surface enters the lower hemisphere at approximately 
1.5 seconds and the collision occurs at approximately 2.0 seconds, at which 
point the material near the collision region begins to compress. We find that 
the conditions at the ignition point are insensitive to the simulation resolution. 

eruption point. We have performed two-dimensional simula- 
tions of various resolutions that begin from a 16 km radius 
burned region offset from the center by 40 km. The history 
of the maximum temperature, Tmax, in the lower hemisphere 
(9 > tt/2), and the density at that same point is shown in Fig- 
ure |5] After the bubble has reached the surface, material - 
burned and unburned - begins to flow over the surface of 
the star towards the opposite pole. This material enters the 
lower hemisphere approximately 1.5 seconds after the igni- 
tion. Some of the low-density material is being shocked as 
it interacts with the stellar surface as it is moving around the 
star, giving a temperature near 10^ K. The collision of mate- 
rial at the lower pole occurs at approximately 2 seconds and 
the density of the hottest region steadily increases thereafter 
until the expected ignition of the detonation at T > 3 x 10^ K 
and p > 10^ g cm"^. See section |5] for a more extensive dis- 
cussion of the fluid motions in the simulation. We find that 
the temperature and density reached at the possible ignition 
point of the detonation are insensitive to the simulation res- 
olution. This gives confidence in the robustness of the GCD 
scenario as a whole, in that the bulk motion of the surface flow 
seems to be insensitive to resolution. Sensitivity to initial con- 
ditions and assumed symmetry have not been addressed here, 
and will be the subject of future work. 

The overall amount of material burned, and therefore the 
amount of energy added to the star, is important for determi- 
nation of the interior state (notably density) when the deto- 
nation wave sweeps through. This will set the amount and 
distribution of ^^Ni and intermediate mass elements, and is 
therefore extremely important for the predictive power of 
our simulations. Shown in Figure |6] is the amount of mass 
burned by the rising flame and the mass of the star, which 
has p > 5.5 X 10^ (see next section). As seen in Figure]?] it 




time (s) 



Fig. 6. — Burned mass (as a fraction of the star) and stellar mass with 
p > 5.5 X lO' g cm"' for simulations shown in figure Js] There is a clear 
dependence on resolution, with finer resolutions generally leading to more 
rapid evolution. 

appears that the bubble rises somewhat faster at higher reso- 
lutions, possibly due to the decreased numerical viscosity or 
to faster burning during the resolved R-T stage. Both the total 
burned mass and the amount of high density material show 
significant dependence on resolution, but it appears that con- 
vergence in the overall quantities may be just within reach. 
The total burned mass at 2 km and 1 km resolution is quite 
consistent by the time the detonation might occur, and the dif- 
ference between the amount of high density material is also 
fairly modest considering the factor of 2 difference in resolu- 
tion and the expected offset in time due to the lower numerical 
viscosity. 

But not all the news is good. Performing a full suite of sim- 
ulations at 1 km was not undertaken for this study and will 
be prohibitive in three dimensions, but knowing how far away 
convergence is lends great interpretive power to our lower res- 
olution simulations. There is also a caution to be raised. Fig- 
ure ]6] does not include results at 0.5 km, although the begin- 
ning of the curve is obviously available, because much more 
mass is burned in this case. This occurs because while the 
dominant trailing vortices follow the main bubble out of the 
star at all lower resolutions, at 0.5 km the branch feature vis- 
ible in Figure ]4] at a radius of 600 km does not. This piece 
of flame, which has now become a ring due to the axisymme- 
try, stays at high density and continues to burn a significant 
portion of the star It is unfortunately not possible for us to 
judge whether this is realistic or largely an artifact of the very 
different nature of vorticity conservation in two dimensions. 
We do note that simulations with a smaller bubble (2 km) and 
slightly larger offset (60 km) do not show this anomaly, and 
initially proceed much as those at lower resolution. This effect 
deserves close scrutiny as more simulations are performed, es- 
pecially in 3-d. Also the shed vortices would likely not have 
an adverse impact if we were simulating a centrally ignited 
deflagration that consumed most of the WD on the way to the 
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Fig. 7. — Maximum temperature Tmax in the lower hemispliere and the 
density at the same point. Shown are results from simulations at 4 km reso- 
lution with an ignited bubble of 16 km radius placed at several offsets from 
the center of the star, 20, 40, 50, 80, and 100 km. Material flowing over 
the surface enters the lower hemisphere at approximately 1.5 seconds and the 
collision occurs at approximately 2.0 seconds, varying some with offset, at 
which point the material near the collision point begins to compress. 

surface. 

5. THE DISTRIBUTIONS OF OUTCOMES PRIOR TO POSSIBLE 
DETONATION 

The uncertainty in the location of the initial flamelet, dis- 
cussed in section |4] leads us to consider ignition of a flame 
bubble at several offsets, roff, from the center of the star We 
find that both the time at which the detonation conditions are 
reached at the point opposite bubble eruption and the expan- 
sion of the star up to this time are correlated with roff. Since 
the expansion of the star is directly related to the density of 
the material through which the detonation will propagate, the 
observed result will be a variation in the ^^Ni mass ejected, 
and that of intermediate elements. 

The maximum temperature, Imax, in the lower hemisphere 
and the density at the same point are shown in Figure [T] for 
a range of offsets from 20 to 100 km. The rise in tempera- 
ture near f = 1 .5 seconds is when the material flowing over 
the surface of the star passes the equator Collision of the 
surface flow at the lower pole occurs at a variety of times be- 
tween 1.7 and 2.4 seconds, where 7n,ax rises to several lO' 
K, and the density steadily increases. Conservative detona- 
tion conditions require T > 3 x 10*^ K and p > 10^ g cm"^ 
dRopke et al.ll2006bt iNiemever & Wooslevl[T997l) . These are 
met at 2.02, 2.05, 2.19 and 2.31 seconds for 100, 80, 50 and 
40 km respectively. Several other properties of the star at the 
time when the ignition is expected to occur are listed in Ta- 
ble |2] particularly the total energy released up to that point, 
in addition to the total burned mass. Note that the total bind- 
ing energy is 4.95 x 10^' erg, so that none of cases here come 
close to unbinding the star; that is expected to occur during 
the detonation phase. Having only flame burning, our models 
continue after the detonation point and show that the detona- 
tion conditions appear to be robust and long-lived, especially 




9(} 



time (s) 



Fig. 8. — Radius and polar angle, 9, location of Tlnax point for simulation 
with rbub = 16 km and roff = 40 km. The region above the surface (r ~ 2 X 
10** cm) is heated as material collides. Material pushed toward the star from 
the collision point then interacts with the stellar surface to form a hot, dense 
region that penetrates into the outer layers of the star 

at the larger offsets. Our 20 km simulation expands the star 
much more by the time of the collision, due to more burning 
occurring in the interior, and may not reach conditions suf- 
ficient for detonation. This certainly indicates that the GCD 
mechanism is likely to fail for ignition points very close to the 
center. 

The radius and the polar angle {0) of the maximum temper- 
ature point are shown in Figure[8j demonstrating the evolution 
of conditions that lead to the detonation ignition. The initial 
surface of the star is at r = 2 x 10^ cm. Thus we see that 
the temperature maximum moves around the surface between 
1.5 and 2.0 seconds, then shifts to the pole at the collision. 
Note that this jump is not material motion. At approximately 
2.2 seconds, material confined between the collision point and 
the star becomes the hottest, and the hot spot moves steadily 
inward from 2x10** cm. During this time the hot spot is not 
always precisely on the axis. Eventually the compression sub- 
sides and this hot spot dissipates. While the hot spot moves 
into the star with a speed of approximately 10** cm s"', ma- 
terial ahead of it (closer to the star) is nearly at rest, and that 
behind it is moving in at just above 10"^ cm s"', so that the hot 
spot occurs in the accumulation. 

A detailed view of the flow near the collision region is 
shown in Figure |9] where we see that a stagnation point is 
formed in the colliding unburned material. From this, mate- 
rial is projected out along the axis and toward the stellar sur- 
face, compressing the surface layers of the star toward ignition 
conditions. The phrase "gravitational confinement" does not 
convey the full impression of what is occurring. The detonat- 
ing material appears to be inertially confined by flow origi- 
nating at the collision point, although the amount of compres- 
sion occurring likely reflects both the strongly gravitationally 
stratified WD surface and a certain amount of assistance from 
gravity, such that both high gravity and and a flow with sig- 
nificant inertia are required to reach such high temperatures 
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TABLE 2 

Properties at Detonation Ignition 



''bub '"off resolution M^,^^^ mass at high density max density energy release 

(km) (km) (km) (s) (% of star) (Mq) (lO'^gcm-^) (10^" g^g-, 



Resolution Study 



16 


40 


8 


2.35 


1.97 


1.04 


7.6 


0.329 


16 


40 


4 


2.31 


2.33 


1.01 


6.6 


0.388 


16 


40 


2 


2.37 


2.90 


0.926 


5.0 


0.478 


16 


40 


1 


2.35'' 


2.88 


0.892 


4.5 


0.517 



Offset Study 



16 


100 


4 


2.02 


1.13 


1.13 


12 


0.193 


16 


80 


4 


2.05 


1.36 


1.12 


11 


0.226 


16 


50 


4 


2.19 


2.07 


1.05 


8.0 


0.347 


16 


40 


4 


2.31 


2.33 


1.01 


6.6 


0.388 


16 


20 


4 


2.70'^ 


6.57 


0.473 


1.6 


1.15 



Note. — All values are evaluated at the time indicated, t^^f. 
" fdet is defined as the first time at which p > 10^ g cm"^^ at the point of Tmux > 3 X lO' K.'' In this case we 
have neglected the early, short-lived, fluctuation at ? ~ 2.15 s"^ Our conservative detonation criteria are not 
reached, listed are values for the peak density of the TiTiax point. 



and densities. The collision itself arises because the mate- 
rial is gravitationally bound, however, it is the kinetic motion 
imparted to the material by the expanded bubble at the break- 
out point that eventually leads to the (gravitationally assisted) 
confinement. 

In the GCD scenario, because so little material is burned 
during the deflagration phase, the amount of ^''Ni produced in 
the supernova is determined by the density distribution during 
the detonation phase. In lieu of simulating the propagation of 
the detonation, which will be performed in future work, we 
have measured the mass of material above 5.5 x 10^ g cm"-^. 
This limit is obtained from the density at which material in 
the W7 model (Nomoto et al. 1984) burned to only 50% ^^Ni. 
This is obviously only a rough estimate, but is good enough 
for measuring the trend with offset distance that we are in- 
terested in here. The bottom panel of Figure [TO] shows how 
this possible ^^Ni mass decreases as the star expands during 
the deflagration phase. The curves are marked at the expected 
launch time of the detonation, where the temperature and den- 
sity first exceed 3 x 10'^ K and 10^ g cm"-' together. 

We find that the amount of ^^Ni expected in the ejecta is 
correlated with the offset of the initial (small) ignition region. 
Larger offsets can produce more ^^Ni for two reasons: ( 1 ) less 
energy is released in the deflagration phase, and therefore the 
star has expanded less when the detonation occurs, and (2) the 
detonation conditions happen sooner so that the star has had 
less time to expand. It does appear that the first of these is the 
dominant effect. The top panel of Figure [TO] shows the mass 
burned as a fraction of the star with time. Larger offsets burn 
less of the star during the bubble rise and breakout, leading to 
less expansion of the star. 

The ^^Ni mass estimates we have found here are fairly high, 
but as seen in section 14.31 this resolution (4 km) appears to 
somewhat underestimate the burned mass and overestimate 
the possible ^^Ni at the detonation time. We are not claiming 
to have performed an absolute calculation of the ^^Ni mass for 
a given ignition point offset; we have instead demonstrated 
a trend that appears to be robust with respect to the physi- 
cal processes that are occurring. We hope that in the future, 
with higher resolution such that self-regulation of the burning 
is strong enough that we can constrain the R-T phase better, 
we may be able to construct a fully predictive model. There 



are, however several steps that should be taken in the mean 
time, including three-dimensional studies that are underway 
( Jordan et al. 2007), and studies of flame bubble response to 
the strong convection expected to be present in the WD core 
when the ignition occurs. The current level of calculation is, 
however sufficient for measuring trends such as how things 
might change with the relative C/O fraction in the interior of 
the WD. 

6. CONCLUSIONS 

We have shown that in the GCD picture of a delayed deto- 
nation of a WD near the Chandrasekhar mass, the properties 
of the WD at detonation, notable the density distribution, are 
systematically correlated with the offset of the ignition point 
of the deflagration. Assuming that the detonation phase pro- 
ceeds as in previous simulations, this will cause a variation in 
the ^^Ni mass ejected in the supernova. The position of the 
ignition point within the inner few 100 km of the WD is ex- 
pected to be stochastically determined by the turbulent flow 
in this region. GCD thus provides a possible explanation for 
the variety of ^^Ni masses seen in Type la Supernovae. 

We find that the conditions (temperature and density 
reached) at the candidate launch point of the detonation are 
insensitive to the resolution of the simulation for resolutions 
studied here (< 8 km). This is a good mark for the robustness 
of the GCD mechanism, but more work is needed, especially 
related to the possibility of vortex shedding early in the bub- 
ble rise and the strong convection that should be present in the 
core at the time of ignition. We have indications of numeri- 
cal convergence in both the total burned mass and the mass 
of dense material, and therefore the predicted ^^Ni mass pro- 
duced by a given ignition offset. But caution is advisable: the 
mass burned during the highly Rayleigh-Taylor (buoyancy- 
driven) unstable rise of the burned region through the star is 
seen to vary with resolution, generally progressing faster with 
higher resolution, even though convergence in the final value 
appears to have been reached. Also, converged results (in the 
extremely limited sense indicated here) appear to require 2 
km or possibly 1 km resolution, which is prohibitive in three 
dimensions. Even here, our parameter study has been per- 
formed at 4 km resolution for efficiency. Thus we are able to 
predict trends in the ^^Ni mass, but not the actual value ejected 
for a given offset. 
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Fig. 9. — Detail of flow near the collision and detonation region for the 
Toff — 40 km case, from top to bottom at ; = 2.07, 2. 19, and 2.32 seconds. 
Temperature is shown in color and contours are shown at p = lO' g cm"'' 
(blue) and at the edge of the burned material ((^ = 0. 1, red). Velocity vectors 
smaller that 10** cm s"' are not shown. A stagnation point is formed above 
the surface of the star from which material is projected out along the axis and 
compressed against the surface of the star, where the detonation is expected 
to occur 




1 1.5 
time (s) 

Fig. 10. — Burned mass (as a fraction of the star) and stellar mass with 
p > 5.5 X 10^ g cm"' for same simulations as in Figure|2] The time at which 
a detonation is expected to be launched is marked with a X for each case. 
Larger offsets are expected to produce more ^*Ni in the ejected material. 



Our method for following the nuclear energy release, in- 
cluding neutronization, with an ADR flame model was de- 
scribed in detail. This method reproduces the energy release 
and hydrodynamic characteristics of the nuclear burning by 
following a limited number of parameters coupled to an arti- 
ficially thickened flame front. We have demonstrated that the 
energy release adds a minimal amount of unwanted acoustic 
noise (RMS velocity < few x 10""*) to the simulation, largely 
removing this source of unrealistic seeds for the instabilities 
in the rising flame surface. 
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